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Abstract 

A new approach to Dynamical Monte Carlo Methods is introduced to simulate markovian pro- 
cesses. We apply this approach to formulate and study an epidemic Generalized SIRS model. The 
results are in excellent agreement with the forth order Runge-Kutta Method in a region of deter- 
ministic solution. We also demonstrate that purely local interactions reproduce a poissonian-like 
process at mesoscopic level. The simulations for this case are checked self-consistently using a 
stochastic version of the Euler Method. 
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I- Introduction - Monte Carlo (MC) methods have been used mainly to equilibrium 
systems]!]], and they have broad applications, since simple systems like hard spheres[0] up 
to complex systems like proteins0, ||. Good reviews in applications of MC methods to 
statistical physics can be seen in the references [pi [5|. In the last decades the development 
of techniques dealing with non-equilibrium systems has been increased 0, specially those 
that concern with stochastic processes. Several attempts were done [[J -JIT] to simulate real 
time processes with this method. Some success was achieved within the scope of poissonian 
processes [|J that has been only recently properly formalized by Fichtorn and Weinberg |§ . 
Another important approach from a theoretical point of view is the waiting (or residence) 
time distribution used by Prados et al.P], whose application is limited to simple systems, 



like Ising models. Some improvement in the real time calculation was presented by Cao||10 
but in a particular and non rigorous way. In this letter we surmount this problem using 
directly the Master Equation, ignoring thus what type of distribution we are dealing. In 
this way, we also avoid the direct waiting (fine-grained) time distribution calculation; this 
is substituted by the calculation of interevent (coarse-grained) times. In our approach, the 
time is a dependent stochastic variable whose distribution is constructed from the Master 
Equation with appropriate transition probabilities. This gives the hierarchy of the process. 
The approach is developed for a class of markovian processes with no simultaneous events in 
the smallest scale considered. Thus, it is for a restricted markovian, but more general than 
poissonian processes. This method has already been applied Jl4| to an extensive study of 



the epidemic Susceptible- Infected-Recovered-Susceptible (SIRS) systems (to details of these 
epidemic systems see and references therein). Here, we apply this new approach to 
formulate an epidemic Generalized SIRS (GSIRS) model, and study two particular cases of 
it. 

II-The Method - For discrete systems, the markovian Master Equation is given by: 

^ = X>^, - X>*p 4 , a) 

3 3 

where Pj is the probability to find the system at the state % at the time t, and w^j is the 
transition probability per unity of time. Considering the probability of transition from 
i to j, we may write w^j = -^-[12], where Tj is a time constant (lifetime) characteristic of 
the state i. 

We now start by choosing a convenient physical extensive microscopic quantity A, that 



is time independent for each state i. The mean value for this quantity at the time t is given 
by: 

A(t) = (A) = J2m)A. (2) 

i 

This equation represents a continuous physical macroscopic quantity A(t). We can differ- 
entiate both sides of the equation above, with respect to t. After that, using (1), and by 
defining AAij = A4 — Aj, we get 

dA(t) 



dt 



-^A-l,- (3) 
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Consider now the nearest-neighbor states j of a given state i; if we measure the "distance" 
between the states, say by the quantity |AAjj|, such that the non-null minimum value is 
|AAjj| = a, we may approach the equation(3) by: 

= Wj^iPjaSij, (4) 

<ij> 

where the symbol < ij > denotes a nearest-neighbour pair of states, and 5^ = AA^/I AA^j |. 
Now we consider another physical quantity A^ that is a source for the quantity A. Thus, we 
can rewrite (4) as: 



dA(t) 



E^M-E^ p A' (5) 



dt 
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where Tj =< Wj_>i >, are the transition probabilities per unity of time averaged over the 
ensemble of the nearest-neighbour states i of j at some time t, i.e., the mesoscopic rates. 
Here, ensemble means a set of configurations accessible at a some finite (small) time around 
a time t; in this sense we are using a time dependent ergodicity idea [5], and so generally 
the systems are non ergodic in non equilibrium states. The superscripts " + " and " — " 
mean respectively the contributions to increasing and to decreasing the quantity A(t). In 
the particular case that r+ = r + and rj = r~ are constants (or only function of the time) 
we have: 

^A = r + A \ _ r - A ( 6 ) 
dt 

what is the analogous to the kinetic equation for the first order chemical reaction A) ^ A, 
being A^ and A the respective concentrations of the chemical elements A) and A. The 
equilibrium can be reached by imposing the balance at macroscopic (or mesoscopic) level: 



r + A^ = r~A. This follows immediately if we require the detailed balance, but it is not 
necessary at all [p~3|j . 

We can write the equation (|j) in an approximated form of a discrete integral 

n 

A(t) - A(t Q ) ~ Yl w^iP 3 {t k )a5 i3 M k . (7) 

k=0 <ij> 

Let now be the set of possible Wj^i represented by Vt = {w^i}, being the states i and 
j occurring around a given instant t, and u>™ ax = supP t . The phase space may be divided 
into N parts, in such way that each part may contain only one element of the system. Thus, 
each element of time in the equation (7) may be represented by 

At k = — ] —. (8) 

We can do the approach to the equation A(t) considering n = EN, with £ sweeps over the 
discretized space; in the limit of N —>■ oo we have the exact solution of the equation (4) for 
a given initial condition. 

Monte Carlo Approach -With the considerations above the equation (0) may be 
written in the form: 

A(t) - A(t ) = f; £ fe) (A) P 3 {t k )a5 ir (9) 

fc=0 <ij> \ t k / \ / 

We can create a hierarchical process choosing the probabilities of transition 

W t k 

that reproduce the correct frequencies of events at each time t k to solve (9). This hierarchy 
have subtle differences with an earlier hierarchy introduced by Fichtorn et al0: first in that 
work (mesoscopic) rates were required, while here we primarily use transition probability per 
unity of time. Second, they used a global maximum to the rates, while here we use a more 
local maximum; in recent work||ll|| this was done without a rigorous proof, based only in 



the detailed balance principle applied to a specific case. To carry out the MC procedure, 
an element is selected randomly with a probability ■h, and thus a transition is tried with 
probability given by (10). The space is swept £ times, with the increment of time in each 
MC step (one MC step here, means a single try to change the state of one element of the 
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system) given by (j8|) up to reach a time t. Starting from the same initial conditions for 
the physical quantities, the process may be repeated, and we can get the average quantity 
A(t) at each instant t. We must emphasize that the probabilities Pj are generated by this 
process. As a given state is chosen with its correct probability in a given time, an ideal MC 
procedure leads to 

tN / 1 \ 

A(t) - A(t ) = £(<^\ - ( r ~ A W ( ^iv ) ■ (") 

where the averages are taken over the ensemble of the states jk at each instant t}~. This is 
just an approach to the integration result of the equation (5). 

We need to observe some important points: first, generally different runs give different 
time tk results at the same MC step k, and the sample averages may be done by linear 
interpolating or extrapolating the data set, in each MC realization, to do them at the same 
point of the time. Second, in one complete sweep around a time tk, the value w™ ax must be 
approximately constant in order do not change the hierarchy and so the result. Third, as the 
configurations do not change drastically in few steps, the microscopic transitions reproduce 
the mesoscopic result. 

Another approach consists in estimating the interevent times by the following rule 

f k a 

At e k = -^f-, (12) 

where rj = r+ and A^ = Aj , or, r| = r~ and A? = Aj k depending on, respectively, 
if the outcome of the experiment increase or decrease the quantity A. The quantity /* is 
an arbitrary e-event dependent factor that must obey the relationship fe = 1j f° r each 

e 

time tk- We emphasize that the time given by (12) represents the average waiting time to 
transitions from a given state jk to any neighbor state i; if the microscopic state remains 
unchanged, the time does not evolve. It can be shown that this procedure leads to the same 
result as using (H) at each MC step observing that 

EE (=$(£) a* < 13 » 

As r e j k A e j k = aY2wj k ^i, using the equation (12) and the normalization condition to in 

i 

(PD , we obtain the expression (|J). In particular, if we choose /* = 0, for most events e, 
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except some e = s, we have /* = 1, so, with this condition, the interevent time has the 
meaning of the waiting time between type-s events. Based on this and in the fact that at 
the equilibrium the relative frequencies of occurrence of events are all equal, we may define 
jk = n k e /Mki where n k e is the number of e— events, and A4 = Yl e n e * s the total number of 
events, in a time interval (arbitrary) near to some time t^. 

III-GSIRS model - Based on (|5|), we formulated the GSIRS model through the 
following set of differential equations and inter-classes rates: 



dS 



dt 



E^ -n'w •»> -i 1 '^- ( 14 ) 



f = E r ^i p A - E ( 15 ) 

3 3 

^ = E ^bPA - E rl^PM' ( 16 ) 

3 3 

where S, I, and R are the populational classes, respectively, of the number of individuals in 
the susceptible, infective and recovered classes. Being the mesoscopic rates r J s ^j , r\^ R and 
r J R ^ s , for each state j, respectively, from S R and R — > S. Note that we meant 

that, for example, if A = I, then A* — S in the equation (|^). The conservation law with 
the total number of individuals iV = S(t) + I(t) + R(t) is satisfied. In particular, a model 
commonly usedO, O give Wr->s = m, ws~,i = T -^S^' 1 ! + A [1 — (1 —po) n ], and w^r = q 
to the transition probabilities per unity of time. We must observe that the mesoscopic rates 
are resulting from local ( "instantaneous" ) averages of the respective transition probabilities 
per unity of time. For practical purposes the individuals are distributed on a square lattice 
of iV = M x M sites. All the individuals at the lattice boundary have their states fixed at 
susceptible state. 

IV-ReSultS and Conclusions - We set the lattice size to M = 200. This 
size was sufficient to get good results compared with the continuum limit when only global 
interactions (A = 0) are considered. The initial condition for the system is set up by 
Iq = 2000 infectives being randomly distributed on the lattice and the remaining sites being 
occupied by So — N — Iq susceptibles, so Ro = 0. 

We consider here two particular cases of the system defined by (|T4| — |T6|) . First, we 
set A = 0, and the other model parameters as q = 0.2, b = 0.8, m = 0.01 and 
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[1 = 2. The non-minimum value, to the differences AS, AI and AR , used in(|T2D is 
a = \AI\ = \AS\ = \AR\ = 1. Figure 1 shows the temporal evolution of I{t). Continu- 
ous lines represent numerical (fourth-order Runge-Kutta) checking solutions for the set of 
differential equations (0 — 0), and open circles correspond to the MC simulations. The 
accuracy of the deterministic solution (Runge-Kutta) was estimated as less than 0.1% (see 
ref.jTT]). Results to the system far from equilibrium showed that the interevent times 
given by (^) have poissonian-like distribution (see inset in figure 1) as expected At the 
equilibrium, the present method leads to converge the distributions of interevent times to 
delta distributions, because the values to the rates and other physical quantities converge 
to constant values. A total of 4 x 10 6 steps, corresponding to 3, 5 x 10 5 configurations, was 
generated by the MC procedure, leading to a total real time of approximately 500 days. 
The total number of configurations used to get the interevent times distribution was about 
8 x 10 4 , what corresponds to approximately 60 days. Second, we set T = and m = 0, i.e., a 
SIR system with purely local variables. The variable n is an integer ranging from n = up 
to 8, since the first and second nearest infected neighbors are indistinguishably considered 
for each susceptible. To this case we use again the expression flT2|), but the rates rs^i are 
obtained by averaging the individual probabilities to the configurations in every successful 
event. This may coast some simulation time. A good optimization for an approximation to 
the exact average is done by drawing randomly susceptibles (1000 here was sufficient) for 
each configuration reached and doing a sample mean with the site transition probabilities 
per unit of time ws^i- It must be observed that this type of average is equivalent to let 
the system advance some small time and take an average over the sample. As the system 
configurations do not change much around some time t^, the small time average corresponds 
to an average in an instantaneous time. To see the self-consistency of the approach, we 
integrate numerically (|14| — ^) given constant (or piecewise constant) time step as in ([?D 
by choosing the maximum local transition probability per unit of time. This maximum is in 
fact actualized at each MC step, when necessary, using a table. When a transition changes a 
state of an individual that changes the maximum, the table is updated. The quantities S, I 
and R are calculated with iterations; the rates are chosen randomly by the MC procedure, 
and thus we use the Euler Method procedure to solve first order differential equations. Ex- 
periments using poissonian distributions^] to obtain the interevent times showed that the 
processes are poissonian-like to all ranges of po, being so, unnecessary the hypothesis of low 
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Po ("weak interaction") as done by Aiello et al|TT|. To illustrate, we show in the Figure 2 
the results to po = 0.8. We compare, also, in figure 2 the iterative method with the MC 
technique described above (restricted markovian method), estimating the interevent time by 
flT2|). The total number of configurations used in the MC procedure was about 4 x 10 4 what 
gives approximately 10 days. The results are in excellent agreement among them. For both 
cases (Figures 1 and 2), results with respect to the MC simulations correspond to an average 
of 20 independent trajectories. The typical MC data errors are in the interval 0.1-1.0%, so 
most of the error bars are smaller than the symbols in the figures. 

We believe that the class of epidemic SIRS models studied here are poissonian-like in 
the mesoscopic scale because of two factors. First, the approach itself implies that no two 
or more events occur in a short scale of time. Second, the mesoscopic rates are slowly 
varying with the time, resembling the independence between events. So, the two conditions 
for a poissonian process were met. We emphasize that low correlations between events are 
not required. It is necessary that the results for independent runs be uncorrelated, so we 
can use the averages obtained for each time t to represent properly the physical quantities 
of the process. To do this we need a local equilibrium hypothesis, what may be at first 
glance restrictive, however we may even reduce the time observation sufficiently such that 
the system does not have time to leave some metastable states. So, we can average it there. 
In the practice of the simulation this is done by increasing the number of observations, i.e., 
the number of time experiments. In forthcoming works we expect to generalize still more 
the method, including up to non-markovian processes. 

The authors gratefully acknowledges funding support from FAPESP Grant n. 00/11635-7 
and 97/03575-0. The authors would also like to thank Drs. F.L.B. da Silva and A. Caliri 
for many stimulating discussions and suggestions. 
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Figure Captions 



FIG. 1. Infected numbers I(t) vs Time. Continuos line: numerical forth-order Runge- 
Kutta solution. Open circles: restricted markovian DMC simulation. Inset: shows the 
behavior of the interevent time At distribution. 

FIG. 2. Infected numbers I(t) vs Time. Continuos line: Iterative stochastic Euler Method 
solution. Squares: restricted markovian DMC simulation. Open circles: poissonian DMC 
simulation. 
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Figure 1 
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Figure 2 
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